Dependencies
library(ggridges) # ggarrange
library(ggpubr)
library(ggdist)
library(ggtext)
library(ggforce)
library(tidyverse) # data manipulation
library(modelr) # data manipulation
library(rstan) #detectCores
library(brms) # core pack
library(parallel) #detectCores
library(tidybayes)
# a flag to let us load models from the disk (F) or re-run the models (T)
# it took about 20 ~ 60 minutes to run the models
REBUILD_FLAG <<- FALSEmetas <- read.csv("../data/metas.csv", stringsAsFactors = F, header = T)
trials_all <- read.csv("../data/trials.csv", stringsAsFactors = F, header = T)
seper <- " o "
trials_all <- trials_all %>%
filter(TrialStatus == "main") %>% #drop practice trial
mutate(Error = TrialAnswer - TrialIntAccuracy,
ErrorMagnitude = abs(Error),
ResponseTime = (TrialVisualizationEnd - TrialVisualizationStart) / 1000,
LogResponseTime = log(ResponseTime)) %>% # log transform
filter(!(ParticipantID == "P7"
& SessionIndex == 0
& TrialStatus == "main"
& TrialIndex == 16)) %>% # filter mistakes
filter(!(ParticipantID == "P9"
& SessionIndex == 2
& TrialStatus == "main"
& TrialIndex == 17)) %>%
filter(!(ParticipantID == "P11"
& SessionIndex == 0
& TrialStatus == "main"
& TrialIndex == 12)) %>%
filter(!(ParticipantID == "P20"
& SessionIndex == 3
& TrialStatus == "main"
& TrialIndex == 5)) %>%
filter(!(ParticipantID == "P24"
& SessionIndex == 0
& TrialStatus == "main"
& TrialIndex == 26)) %>%
filter(ErrorMagnitude < 51) %>%
mutate(
# manual recode levels using orthogonal contrast coding
# the distance within a variable is always 1
AllFactors = paste(FactorStereo, FactorMotion, FactorPerspective, FactorShading, Dim, sep = seper),
Device = recode(FactorStereo, yes = 0.5, no = -0.5),
Motion = recode(FactorMotion, yes = 0.5, no = -0.5),
Projection = recode(FactorPerspective, yes = 0.5, no = -0.5),
# this will ignore shading-1, the commonly-used phong lighting model
Shading = recode(FactorShading, `shading-1` = 0, `shading-2` = 0.5, `flat`= -0.5, `no` = -0.5),
DataModel = recode(TrialDataset, "cifar10" = 0.5, "babi-q1" = -0.5),
Dimensionality = recode(Dim, `2` = -0.5, `3` = 0.5)
)
trials <- trials_all %>%
select(ParticipantID, TrialIntAccuracy, ErrorMagnitude, ResponseTime, Device, Motion, Projection, Shading, Dimensionality, DataModel)
vtrials <- trials_all %>%
select(ParticipantID, ErrorMagnitude, ResponseTime, FactorStereo, FactorMotion, FactorPerspective, FactorShading, Dim, TrialDataset, AllFactors)We first plot all raw data to get a sense of how data distribute. Both error magnitude and response time datum seem not to follow a normal distribution, and they are both non-negative. A skewed normal, lognormal, and gamma distribution could be reasonable. Note that error magnitude contains 0, and response time has a upper bound 60s (e.g., 60.1s, 60.2s due to running time). We will use gamma distribution for both measures. To formally compare which distribution better fits data, we could use the model_diagnostics function defined below.
g_raw_error <- vtrials %>%
ggplot(mapping = aes(x = ErrorMagnitude)) +
facet_wrap(AllFactors ~ .,ncol = 3, scales = "free") +
coord_cartesian(xlim = c(0, 50), ylim = c(0, 30)) +
geom_histogram(binwidth = 2, fill = "lightblue", color = "white") +
ggtitle("raw data - error magnitude (pt)") +
xlab ("error magnitude ") + ylab("count")
ggsave(g_raw_error, filename = "pic/g_raw_error.pdf", width = 10, height = 10, units = "in")
g_raw_errorg_raw_time <- vtrials %>%
ggplot(mapping = aes(x = ResponseTime)) +
facet_wrap(AllFactors ~ .,ncol = 3, scales = "free") +
coord_cartesian(xlim = c(0, 62), ylim = c(0, 30)) +
geom_histogram(binwidth = 2, fill = "lightgray", color = "white") +
ggtitle("raw data - response time (s)") +
xlab ("response time ") + ylab("count")
ggsave(g_raw_time, filename = "pic/g_raw_time.pdf", width = 10, height = 10, units = "in")
g_raw_time
# Models
First we define a helper function to wrap up model diagnostics.
model_diagnostics <- function(m, name = "model diagnostics"){
print(name)
print(summary(m))
# plot(m) # run this line to check convergence when first fitting the model
plot(pp_check(m, type = "dens_overlay", nsamples = 100))
plot(pp_check(m, type = "stat_2d", nsamples = 2000))
plot(pp_check(m, type = "error_binned", nsamples = 100))
print(loo(m))
print(waic(m))
}For error magnitude, we have 0. Therefore, we use a hurdle_gamma distribution to account for this. Based on the experimental design, Dimensionality only interacts with Device and DataModel. The other five variables could interact with each other. We are also concerned about the difference between Dimensionality; we therefore model it as a group-level effect. The two data models have different ranges of accuracy; we also specify it as a group-effect. Last, each trial has different true accuracies, and we specify it as another random intercept to separate the effects. Because of the skewness (we also don’t have evidence of error magnitude being normally distributed), we use a gamma likelihood instead of normal/student_t. Because of the zero values (no error), we use a hurdle gamma likelihood. Other choices like skewed normal are also reasonable, but gamma results in posterior prediction best capturing the data and of a smaller waic. Since we care about mean of error magnitude and the data also does not suggest a big variance in standard deviation across conditions, we did not use a submodel for standard deviation.
if(REBUILD_FLAG){
m_error_hurdle_gamma <- bf(ErrorMagnitude ~
0 + Intercept +
Device*Motion*Shading*Projection*DataModel +
Device*Dimensionality*DataModel +
(1 + Dimensionality*DataModel | ParticipantID) +
(1 | TrialIntAccuracy),
family = hurdle_gamma(link = "log")
)
get_prior(m_error_hurdle_gamma, data = trials)
m_brms_error_hurdle_gamma <<- brm(m_error_hurdle_gamma,
data = trials,
family = hurdle_gamma(link = "log"),
prior = c( # priors" 2 SDs roughly catch 95% data
prior(normal(0, 2.5), class = "b"),
prior(normal(0, 25), class = "b", coef = "Intercept"),
prior(gamma(2, 10), class = "shape"),
prior(student_t(3, 0, 2.5), class = "sd")
),
backend = 'cmdstanr',
chains = N_CHAINS, iter = N_ITER, warmup = N_WARMUP, thin = N_THIN, cores = N_CORES,
control = list(adapt_delta = DELTA, max_treedepth = TREE_DEPTH))
saveRDS(file = "rds/m_brms_error_hurdle_gamma.rds", m_brms_error_hurdle_gamma)
}else{
m_brms_error_hurdle_gamma <<- readRDS("rds/m_brms_error_hurdle_gamma.rds")
}
model_diagnostics(m_brms_error_hurdle_gamma)## [1] "model diagnostics"
## Family: hurdle_gamma
## Links: mu = log; shape = identity; hu = identity
## Formula: ErrorMagnitude ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Dimensionality * DataModel | ParticipantID) + (1 | TrialIntAccuracy)
## Data: trials (Number of observations: 3320)
## Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~ParticipantID (Number of levels: 32)
## Estimate Est.Error l-95% CI
## sd(Intercept) 0.21 0.04 0.14
## sd(Dimensionality) 0.14 0.08 0.01
## sd(DataModel) 0.38 0.07 0.25
## sd(Dimensionality:DataModel) 0.20 0.15 0.01
## cor(Intercept,Dimensionality) 0.15 0.37 -0.58
## cor(Intercept,DataModel) -0.16 0.23 -0.59
## cor(Dimensionality,DataModel) -0.35 0.34 -0.88
## cor(Intercept,Dimensionality:DataModel) -0.08 0.41 -0.81
## cor(Dimensionality,Dimensionality:DataModel) -0.13 0.43 -0.86
## cor(DataModel,Dimensionality:DataModel) 0.10 0.41 -0.72
## u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.30 1.00 2461 3523
## sd(Dimensionality) 0.30 1.00 1341 2692
## sd(DataModel) 0.54 1.00 2566 3284
## sd(Dimensionality:DataModel) 0.54 1.00 1830 2654
## cor(Intercept,Dimensionality) 0.82 1.00 3101 3255
## cor(Intercept,DataModel) 0.32 1.00 1920 2239
## cor(Dimensionality,DataModel) 0.40 1.00 661 1023
## cor(Intercept,Dimensionality:DataModel) 0.72 1.00 3422 3497
## cor(Dimensionality,Dimensionality:DataModel) 0.73 1.00 2625 3253
## cor(DataModel,Dimensionality:DataModel) 0.83 1.00 3607 3287
##
## ~TrialIntAccuracy (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.31 0.04 0.24 0.40 1.00 1935 2886
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 2.02 0.07 1.88 2.16
## Device -0.01 0.07 -0.15 0.12
## Motion -0.14 0.04 -0.21 -0.06
## Shading -0.03 0.03 -0.10 0.04
## Projection -0.01 0.03 -0.06 0.05
## DataModel 0.06 0.10 -0.14 0.25
## Dimensionality 0.10 0.07 -0.04 0.25
## Device:Motion 0.01 0.07 -0.14 0.16
## Device:Shading -0.01 0.07 -0.15 0.12
## Motion:Shading 0.08 0.07 -0.05 0.22
## Device:Projection -0.03 0.06 -0.15 0.08
## Motion:Projection -0.04 0.06 -0.15 0.07
## Shading:Projection 0.06 0.07 -0.08 0.19
## Device:DataModel 0.12 0.14 -0.15 0.38
## Motion:DataModel -0.09 0.07 -0.24 0.05
## Shading:DataModel -0.05 0.07 -0.19 0.09
## Projection:DataModel 0.08 0.06 -0.03 0.20
## Device:Dimensionality 0.13 0.14 -0.13 0.41
## DataModel:Dimensionality 0.13 0.14 -0.16 0.40
## Device:Motion:Shading 0.07 0.14 -0.21 0.34
## Device:Motion:Projection 0.07 0.11 -0.14 0.29
## Device:Shading:Projection 0.22 0.14 -0.04 0.48
## Motion:Shading:Projection 0.04 0.14 -0.23 0.33
## Device:Motion:DataModel 0.06 0.15 -0.23 0.35
## Device:Shading:DataModel 0.02 0.14 -0.25 0.29
## Motion:Shading:DataModel 0.40 0.14 0.12 0.66
## Device:Projection:DataModel 0.14 0.11 -0.07 0.36
## Motion:Projection:DataModel -0.01 0.11 -0.23 0.21
## Shading:Projection:DataModel -0.12 0.14 -0.40 0.15
## Device:DataModel:Dimensionality -0.12 0.27 -0.65 0.40
## Device:Motion:Shading:Projection 0.00 0.28 -0.53 0.55
## Device:Motion:Shading:DataModel 0.24 0.26 -0.27 0.76
## Device:Motion:Projection:DataModel -0.20 0.23 -0.65 0.25
## Device:Shading:Projection:DataModel -0.01 0.27 -0.54 0.53
## Motion:Shading:Projection:DataModel 0.45 0.27 -0.10 0.98
## Device:Motion:Shading:Projection:DataModel -0.10 0.55 -1.19 0.95
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1954 3030
## Device 1.00 3520 3500
## Motion 1.00 3542 3488
## Shading 1.00 3429 3319
## Projection 1.00 3491 3455
## DataModel 1.00 3159 3204
## Dimensionality 1.00 3189 3462
## Device:Motion 1.00 3465 3294
## Device:Shading 1.00 3442 3316
## Motion:Shading 1.00 3552 2685
## Device:Projection 1.00 3244 3246
## Motion:Projection 1.00 3258 3183
## Shading:Projection 1.00 3360 3312
## Device:DataModel 1.00 3770 3464
## Motion:DataModel 1.00 3505 3692
## Shading:DataModel 1.00 3368 3439
## Projection:DataModel 1.00 3441 2814
## Device:Dimensionality 1.00 3593 3604
## DataModel:Dimensionality 1.00 3131 3450
## Device:Motion:Shading 1.00 3595 3334
## Device:Motion:Projection 1.00 3551 3228
## Device:Shading:Projection 1.00 3297 3229
## Motion:Shading:Projection 1.00 3468 3099
## Device:Motion:DataModel 1.00 3667 3435
## Device:Shading:DataModel 1.00 3470 3180
## Motion:Shading:DataModel 1.00 2700 3165
## Device:Projection:DataModel 1.00 3083 2932
## Motion:Projection:DataModel 1.00 3354 3284
## Shading:Projection:DataModel 1.00 3367 3183
## Device:DataModel:Dimensionality 1.00 3883 3283
## Device:Motion:Shading:Projection 1.00 3362 3072
## Device:Motion:Shading:DataModel 1.00 3562 3169
## Device:Motion:Projection:DataModel 1.00 3047 3315
## Device:Shading:Projection:DataModel 1.00 3412 3486
## Motion:Shading:Projection:DataModel 1.00 3524 3156
## Device:Motion:Shading:Projection:DataModel 1.00 3376 2850
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 1.80 0.04 1.72 1.89 1.00 3234 2612
## hu 0.07 0.00 0.06 0.08 1.00 3309 3104
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##
## Computed from 4000 by 3320 log-likelihood matrix
##
## Estimate SE
## elpd_loo -9723.0 49.6
## p_loo 129.7 4.6
## looic 19446.0 99.2
## ------
## Monte Carlo SE of elpd_loo is 0.3.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 3316 99.9% 530
## (0.5, 0.7] (ok) 4 0.1% 351
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Computed from 4000 by 3320 log-likelihood matrix
##
## Estimate SE
## elpd_waic -9721.4 49.5
## p_waic 128.2 4.4
## waic 19442.9 99.0
##
## 32 (1.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
For response time, we have a similar motivation with error magnitude with a few exceptions. First response time is positive. Therefore a gamma distribution is good enough. Second, we know the upper bound of response time is about 60 is cropped. We need to tell the model about this (left censored). Also, data model and dimensionality probably does not affect response time very much; however, different devices and motion levels should be treated differently. So we use Device*Motion as the group-level effects to replace data properties as group-level effects. We note that the distributions of response time could look very different; a sub-model for shape is more appropriate to capture the difference. We added Device and Motion as fixed effects (adding random effects would let the model not converge.)
if(REBUILD_FLAG == TRUE){
trials <- trials %>%
mutate(cen = ifelse(ResponseTime > 60, "left", "none"))
m_time_gamma <- bf(
ResponseTime | cens(cen) ~ 0 + Intercept +
Device*Motion*Shading*Projection*DataModel +
Device*Dimensionality*DataModel +
(1 + Device*Motion | ParticipantID),
shape ~ 0 + Intercept +
Device*Motion,
family = Gamma(link="log")
)
priors_time_gamma <- c(
prior(normal(0, 2.5), class = "b"),
prior(normal(0, 30), class = "b", coef = "Intercept"),
#prior(gamma(1, 20), class = "shape"),
prior(student_t(3, 0, 2.5), class = "sd")
)
m_brms_time_gamma <<- brm(m_time_gamma,
data = trials, family = Gamma(link="log"), prior = priors_time_gamma,
backend = 'cmdstanr',
chains = N_CHAINS, iter = N_ITER, warmup = N_WARMUP, thin = N_THIN, cores = N_CORES,
control = list(adapt_delta = DELTA, max_treedepth = TREE_DEPTH))
saveRDS(file = "rds/m_brms_time_gamma.rds", m_brms_time_gamma)
}else{
m_brms_time_gamma <<- readRDS("rds/m_brms_time_gamma.rds")
}
model_diagnostics(m_brms_time_gamma)## [1] "model diagnostics"
## Family: gamma
## Links: mu = log; shape = log
## Formula: ResponseTime | cens(cen) ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Device * Motion | ParticipantID)
## shape ~ 0 + Intercept + Device * Motion
## Data: trials (Number of observations: 3320)
## Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~ParticipantID (Number of levels: 32)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.64 0.08 0.50 0.83 1.00 1523
## sd(Device) 0.25 0.04 0.18 0.34 1.00 2253
## sd(Motion) 0.29 0.05 0.21 0.39 1.00 2545
## sd(Device:Motion) 0.75 0.11 0.56 1.00 1.00 2499
## cor(Intercept,Device) -0.11 0.19 -0.46 0.26 1.00 3245
## cor(Intercept,Motion) -0.48 0.15 -0.73 -0.14 1.00 2975
## cor(Device,Motion) -0.10 0.20 -0.48 0.29 1.00 2583
## cor(Intercept,Device:Motion) -0.14 0.18 -0.47 0.22 1.00 2874
## cor(Device,Device:Motion) 0.05 0.20 -0.36 0.42 1.00 2015
## cor(Motion,Device:Motion) 0.08 0.19 -0.30 0.44 1.00 1829
## Tail_ESS
## sd(Intercept) 2186
## sd(Device) 2867
## sd(Motion) 3241
## sd(Device:Motion) 2510
## cor(Intercept,Device) 3523
## cor(Intercept,Motion) 3479
## cor(Device,Motion) 3367
## cor(Intercept,Device:Motion) 3201
## cor(Device,Device:Motion) 2950
## cor(Motion,Device:Motion) 2740
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 2.27 0.11 2.05 2.50
## Device 0.37 0.06 0.25 0.50
## Motion 0.59 0.05 0.49 0.70
## Shading 0.01 0.02 -0.04 0.05
## Projection 0.04 0.02 0.00 0.08
## DataModel 0.14 0.07 -0.00 0.28
## Dimensionality 0.12 0.05 0.02 0.21
## Device:Motion 0.37 0.14 0.10 0.65
## Device:Shading 0.09 0.05 -0.00 0.18
## Motion:Shading 0.14 0.05 0.05 0.24
## Device:Projection 0.03 0.04 -0.04 0.11
## Motion:Projection 0.05 0.04 -0.02 0.13
## Shading:Projection -0.07 0.05 -0.17 0.02
## Device:DataModel 0.03 0.14 -0.26 0.31
## Motion:DataModel -0.02 0.13 -0.29 0.23
## Shading:DataModel 0.07 0.05 -0.02 0.17
## Projection:DataModel 0.02 0.04 -0.06 0.09
## Device:Dimensionality -0.06 0.10 -0.26 0.13
## DataModel:Dimensionality 0.05 0.10 -0.13 0.25
## Device:Motion:Shading -0.06 0.09 -0.24 0.12
## Device:Motion:Projection -0.05 0.08 -0.20 0.10
## Device:Shading:Projection 0.07 0.09 -0.11 0.26
## Motion:Shading:Projection 0.10 0.10 -0.09 0.29
## Device:Motion:DataModel 0.39 0.26 -0.11 0.90
## Device:Shading:DataModel -0.02 0.10 -0.20 0.17
## Motion:Shading:DataModel 0.05 0.09 -0.13 0.24
## Device:Projection:DataModel -0.03 0.08 -0.18 0.13
## Motion:Projection:DataModel -0.06 0.08 -0.21 0.10
## Shading:Projection:DataModel 0.01 0.09 -0.17 0.19
## Device:DataModel:Dimensionality -0.06 0.19 -0.43 0.31
## Device:Motion:Shading:Projection 0.24 0.18 -0.13 0.59
## Device:Motion:Shading:DataModel -0.12 0.19 -0.48 0.24
## Device:Motion:Projection:DataModel 0.08 0.15 -0.22 0.38
## Device:Shading:Projection:DataModel -0.05 0.19 -0.42 0.32
## Motion:Shading:Projection:DataModel 0.27 0.19 -0.10 0.63
## Device:Motion:Shading:Projection:DataModel 0.20 0.37 -0.51 0.93
## shape_Intercept 1.30 0.02 1.26 1.35
## shape_Device 0.22 0.05 0.13 0.31
## shape_Motion 0.32 0.05 0.22 0.41
## shape_Device:Motion 0.20 0.10 0.01 0.40
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1361 2099
## Device 1.00 3072 3488
## Motion 1.00 2404 3262
## Shading 1.00 3621 3624
## Projection 1.00 3976 3552
## DataModel 1.00 2930 3508
## Dimensionality 1.00 3571 3331
## Device:Motion 1.00 2564 3077
## Device:Shading 1.00 3496 3657
## Motion:Shading 1.00 3564 3428
## Device:Projection 1.00 3704 3397
## Motion:Projection 1.00 3707 3491
## Shading:Projection 1.00 3915 3561
## Device:DataModel 1.00 3009 3055
## Motion:DataModel 1.00 2601 2912
## Shading:DataModel 1.00 3854 3575
## Projection:DataModel 1.00 3657 3362
## Device:Dimensionality 1.00 3402 3159
## DataModel:Dimensionality 1.00 3372 3664
## Device:Motion:Shading 1.00 3559 3334
## Device:Motion:Projection 1.00 3971 3807
## Device:Shading:Projection 1.00 3918 3685
## Motion:Shading:Projection 1.00 3751 3604
## Device:Motion:DataModel 1.00 2701 3242
## Device:Shading:DataModel 1.00 3919 3614
## Motion:Shading:DataModel 1.00 3616 3661
## Device:Projection:DataModel 1.00 3600 3452
## Motion:Projection:DataModel 1.00 3756 3346
## Shading:Projection:DataModel 1.00 3820 3572
## Device:DataModel:Dimensionality 1.00 3543 3879
## Device:Motion:Shading:Projection 1.00 3723 3436
## Device:Motion:Shading:DataModel 1.00 3533 3542
## Device:Motion:Projection:DataModel 1.00 3781 3634
## Device:Shading:Projection:DataModel 1.00 3898 3767
## Motion:Shading:Projection:DataModel 1.00 3679 3492
## Device:Motion:Shading:Projection:DataModel 1.00 3725 3657
## shape_Intercept 1.00 3592 3524
## shape_Device 1.00 3810 3379
## shape_Motion 1.00 3590 3237
## shape_Device:Motion 1.00 3736 3499
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##
## Computed from 4000 by 3320 log-likelihood matrix
##
## Estimate SE
## elpd_loo -9796.9 74.5
## p_loo 179.8 10.5
## looic 19593.9 149.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 3301 99.4% 555
## (0.5, 0.7] (ok) 15 0.5% 155
## (0.7, 1] (bad) 4 0.1% 37
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Computed from 4000 by 3320 log-likelihood matrix
##
## Estimate SE
## elpd_waic -9795.1 74.3
## p_waic 178.0 10.2
## waic 19590.1 148.7
##
## 68 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
summary(m_brms_time_gamma)## Family: gamma
## Links: mu = log; shape = log
## Formula: ResponseTime | cens(cen) ~ 0 + Intercept + Device * Motion * Shading * Projection * DataModel + Device * Dimensionality * DataModel + (1 + Device * Motion | ParticipantID)
## shape ~ 0 + Intercept + Device * Motion
## Data: trials (Number of observations: 3320)
## Draws: 2 chains, each with iter = 2000; warmup = 0; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~ParticipantID (Number of levels: 32)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.64 0.08 0.50 0.83 1.00 1523
## sd(Device) 0.25 0.04 0.18 0.34 1.00 2253
## sd(Motion) 0.29 0.05 0.21 0.39 1.00 2545
## sd(Device:Motion) 0.75 0.11 0.56 1.00 1.00 2499
## cor(Intercept,Device) -0.11 0.19 -0.46 0.26 1.00 3245
## cor(Intercept,Motion) -0.48 0.15 -0.73 -0.14 1.00 2975
## cor(Device,Motion) -0.10 0.20 -0.48 0.29 1.00 2583
## cor(Intercept,Device:Motion) -0.14 0.18 -0.47 0.22 1.00 2874
## cor(Device,Device:Motion) 0.05 0.20 -0.36 0.42 1.00 2015
## cor(Motion,Device:Motion) 0.08 0.19 -0.30 0.44 1.00 1829
## Tail_ESS
## sd(Intercept) 2186
## sd(Device) 2867
## sd(Motion) 3241
## sd(Device:Motion) 2510
## cor(Intercept,Device) 3523
## cor(Intercept,Motion) 3479
## cor(Device,Motion) 3367
## cor(Intercept,Device:Motion) 3201
## cor(Device,Device:Motion) 2950
## cor(Motion,Device:Motion) 2740
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 2.27 0.11 2.05 2.50
## Device 0.37 0.06 0.25 0.50
## Motion 0.59 0.05 0.49 0.70
## Shading 0.01 0.02 -0.04 0.05
## Projection 0.04 0.02 0.00 0.08
## DataModel 0.14 0.07 -0.00 0.28
## Dimensionality 0.12 0.05 0.02 0.21
## Device:Motion 0.37 0.14 0.10 0.65
## Device:Shading 0.09 0.05 -0.00 0.18
## Motion:Shading 0.14 0.05 0.05 0.24
## Device:Projection 0.03 0.04 -0.04 0.11
## Motion:Projection 0.05 0.04 -0.02 0.13
## Shading:Projection -0.07 0.05 -0.17 0.02
## Device:DataModel 0.03 0.14 -0.26 0.31
## Motion:DataModel -0.02 0.13 -0.29 0.23
## Shading:DataModel 0.07 0.05 -0.02 0.17
## Projection:DataModel 0.02 0.04 -0.06 0.09
## Device:Dimensionality -0.06 0.10 -0.26 0.13
## DataModel:Dimensionality 0.05 0.10 -0.13 0.25
## Device:Motion:Shading -0.06 0.09 -0.24 0.12
## Device:Motion:Projection -0.05 0.08 -0.20 0.10
## Device:Shading:Projection 0.07 0.09 -0.11 0.26
## Motion:Shading:Projection 0.10 0.10 -0.09 0.29
## Device:Motion:DataModel 0.39 0.26 -0.11 0.90
## Device:Shading:DataModel -0.02 0.10 -0.20 0.17
## Motion:Shading:DataModel 0.05 0.09 -0.13 0.24
## Device:Projection:DataModel -0.03 0.08 -0.18 0.13
## Motion:Projection:DataModel -0.06 0.08 -0.21 0.10
## Shading:Projection:DataModel 0.01 0.09 -0.17 0.19
## Device:DataModel:Dimensionality -0.06 0.19 -0.43 0.31
## Device:Motion:Shading:Projection 0.24 0.18 -0.13 0.59
## Device:Motion:Shading:DataModel -0.12 0.19 -0.48 0.24
## Device:Motion:Projection:DataModel 0.08 0.15 -0.22 0.38
## Device:Shading:Projection:DataModel -0.05 0.19 -0.42 0.32
## Motion:Shading:Projection:DataModel 0.27 0.19 -0.10 0.63
## Device:Motion:Shading:Projection:DataModel 0.20 0.37 -0.51 0.93
## shape_Intercept 1.30 0.02 1.26 1.35
## shape_Device 0.22 0.05 0.13 0.31
## shape_Motion 0.32 0.05 0.22 0.41
## shape_Device:Motion 0.20 0.10 0.01 0.40
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1361 2099
## Device 1.00 3072 3488
## Motion 1.00 2404 3262
## Shading 1.00 3621 3624
## Projection 1.00 3976 3552
## DataModel 1.00 2930 3508
## Dimensionality 1.00 3571 3331
## Device:Motion 1.00 2564 3077
## Device:Shading 1.00 3496 3657
## Motion:Shading 1.00 3564 3428
## Device:Projection 1.00 3704 3397
## Motion:Projection 1.00 3707 3491
## Shading:Projection 1.00 3915 3561
## Device:DataModel 1.00 3009 3055
## Motion:DataModel 1.00 2601 2912
## Shading:DataModel 1.00 3854 3575
## Projection:DataModel 1.00 3657 3362
## Device:Dimensionality 1.00 3402 3159
## DataModel:Dimensionality 1.00 3372 3664
## Device:Motion:Shading 1.00 3559 3334
## Device:Motion:Projection 1.00 3971 3807
## Device:Shading:Projection 1.00 3918 3685
## Motion:Shading:Projection 1.00 3751 3604
## Device:Motion:DataModel 1.00 2701 3242
## Device:Shading:DataModel 1.00 3919 3614
## Motion:Shading:DataModel 1.00 3616 3661
## Device:Projection:DataModel 1.00 3600 3452
## Motion:Projection:DataModel 1.00 3756 3346
## Shading:Projection:DataModel 1.00 3820 3572
## Device:DataModel:Dimensionality 1.00 3543 3879
## Device:Motion:Shading:Projection 1.00 3723 3436
## Device:Motion:Shading:DataModel 1.00 3533 3542
## Device:Motion:Projection:DataModel 1.00 3781 3634
## Device:Shading:Projection:DataModel 1.00 3898 3767
## Motion:Shading:Projection:DataModel 1.00 3679 3492
## Device:Motion:Shading:Projection:DataModel 1.00 3725 3657
## shape_Intercept 1.00 3592 3524
## shape_Device 1.00 3810 3379
## shape_Motion 1.00 3590 3237
## shape_Device:Motion 1.00 3736 3499
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Now we have the models. Before we head towards guiding questions, we first a helper function for plotting the results.
draw_coefs <- function(m,
target_coefs,
title = "",
fill_color = "white" ,
breaks_x = seq(0, 1, by = 0.5)
){
post_samples <- posterior_samples(m, pars = "b_")
b_CIs <- post_samples %>%
select(target_coefs) %>%
gather(key = "coefs", value = "value") %>%
group_by(coefs) %>%
median_qi() %>%
ungroup() %>%
mutate(lower95 = exp(.lower), upper95 = exp(.upper), median = exp(value))%>%
select(coefs, median, lower95, upper95) %>%
mutate(coefs = factor(coefs, levels = rev(target_coefs)),
s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
CIs = paste0(round(median, 3), ' [', s_lower95, ',', s_upper95, ']'))
b_samples <- post_samples %>%
select(target_coefs) %>%
gather(key = "coefs", value = "value") %>%
ungroup()
pr <- b_samples %>%
mutate(`Pr > 0` = if_else(value >= 0 , 1, 0)) %>%
group_by(coefs) %>%
summarise(`Pr > 0` = sum(`Pr > 0`) / nrow(post_samples)) %>%
mutate(`Pr > 0` = round(`Pr > 0` * 100, 2))
g <- b_samples %>%
mutate(value = exp(value)) %>%
ggplot(aes(x = value, y = coefs)) +
geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.1, y = coefs, label = CIs),
size = 2,
fill = NA, label.color = NA, # remove background and outline
label.padding = grid::unit(rep(0, 4), "pt")) +
geom_richtext(data = pr, mapping = aes(x = tail(breaks_x, n = 1) * 1.25, y = coefs, label = `Pr > 0`),
size = 2,
fill = NA, label.color = NA, # remove background and outline
label.padding = grid::unit(rep(0, 4), "pt")) +
geom_density_ridges(
fill = fill_color, color = fill_color,
rel_min_height = 0.0001,
alpha = 0.5,
scale = .9,
size = 0.5
) +
geom_vline(xintercept = 1, linetype = "dashed", size = LINE_SIZE) +
scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .5))) +
coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
ggtitle(title) + xlab("exp(beta)") + ylab("")
return(g)
}To answer the first guiding question, we look at the model coefficients of each experimental variable from fixed effects. The coefficients from fixed effects tell us the overall weight of each variable when conditioning on just that variable for an average; note that the sign and direction of an effect depends on coding.
target_coefs_rq1 <- c("b_Device", "b_Motion", "b_Projection", "b_Shading", "b_Dimensionality", "b_DataModel")
g1_rq1 <- draw_coefs(m_brms_error_hurdle_gamma,
target_coefs_rq1,
"Error magnitude",
ERROR_COLOR,
breaks_x = seq(0.75, 1.5, by = 0.25))
g2_rq1 <- draw_coefs(m_brms_time_gamma,
target_coefs_rq1,
"Response time",
TIME_COLOR,
breaks_x = seq(0.85, 2.2, by = 0.45))
g_combined_rq1 <- ggarrange(g1_rq1, g2_rq1, ncol = 2)
ggsave(g_combined_rq1, filename = "pic/RQ1-mean-effects.pdf", width = 8, height = 2.25, units = "in")
g_combined_rq1Similar to our investigation for the first guiding question, we use the model coefficients to investigate the interactions effects between variables. We only considered two-way interaction effects for simplicity; more complicated interaction effects were implied by the analyses below in GQ4.
target_coefs_rq2 <- c("b_Device:Motion",
"b_Device:Shading",
"b_Device:Projection",
"b_Device:Dimensionality",
"b_Device:DataModel",
"b_Motion:Shading",
"b_Motion:Projection",
"b_Motion:DataModel",
"b_Shading:Projection",
"b_Shading:DataModel",
"b_Projection:DataModel",
"b_DataModel:Dimensionality")
g1_rq2 <- draw_coefs(m_brms_error_hurdle_gamma,
target_coefs_rq2,
"Error magnitude",
ERROR_COLOR,
breaks_x = seq(0.5, 2.1, by = 0.5))
g2_rq2 <- draw_coefs(m_brms_time_gamma,
target_coefs_rq2,
"Response time",
TIME_COLOR,
breaks_x = seq(0.5, 2.1, by = 0.5))
g_combined_rq2 <- ggarrange(g1_rq2, g2_rq2, ncol = 2)
ggsave(g_combined_rq2, filename = "pic/RQ2-interaction-effects.pdf", width = 8.5, height = 3.75, units = "in")
g_combined_rq2To compare different levels of manipulations, we draw posterior samples from the model to get an estimate of mean for each measure. These posterior distributions show how an average participant does when conditioning on that level of manipulation.
target_order <- c("Device", "Motion", "Projection", "Shading", "Dimensionality", "DataModel")
m <- m_brms_error_hurdle_gamma
draw_post <- function( m,
title = "",
fill_color = "white" ,
palette_in,
adjust_in,
breaks_x = seq(0, 1, by = 0.5)
){
post_samples <- trials %>%
# these use theoritical centers, not practical centers...
# data_grid(Device = c(-0.5, 0.5), Motion = 0, Projection = 0, Shading = 0, Dimensionality = 0, DataModel = 0) %>%
# add_row(trials %>%data_grid(Device = 0, Motion = c(-0.5, 0.5), Projection = 0, Shading = 0, Dimensionality = 0, DataModel = 0)) %>%
# add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = c(-0.5, 0.5), Shading = 0, Dimensionality = 0, DataModel = 0)) %>%
# add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = c(-0.5, 0, 0.5), Dimensionality = 0, DataModel = 0)) %>%
# add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = 0, Dimensionality = c(-0.5, 0.5), DataModel = 0)) %>%
# add_row(trials %>%data_grid(Device = 0, Motion = 0, Projection = 0, Shading = 0, Dimensionality = 0, DataModel = c(-0.5, 0.5))) %>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
distinct() %>%
#View(post_samples)
add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = F) %>%
select(-c(".chain", ".iteration"))%>%
ungroup()%>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
gather(key = "levels_names", value = "levels_values", Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
mutate(levels_names = factor(levels_names, levels = rev(target_order)), levels_values = factor(levels_values))
b_CIs <- post_samples %>%
group_by(levels_names, levels_values) %>%
median_qi() %>%
ungroup() %>%
mutate(lower95 = .lower, upper95 = .upper, median = sample_value)%>%
mutate(levels_values = factor(levels_values),
s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
CIs = paste0(round(median, 3), ' [', s_lower95, ',', s_upper95, ']'))
g <- post_samples %>%
ggplot(aes(x = sample_value, y = levels_names)) +
geom_density_ridges(
aes(fill = NA, color = levels_values),
rel_min_height = 0.0005,
alpha = 0.5,
scale = .8,
size = 0.4
) +
geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.2,
y = levels_names,
color = levels_values,
label = CIs),
label.color = NA,
size = 2,
fill = NA,
label.padding = grid::unit(rep(0, 4), "pt")) +
#scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .5))) +
scale_fill_manual(values = palette_in) +
scale_color_manual(values = palette_in) +
#coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
ggtitle(title) + xlab("") + ylab("")
return(g)
}
g1_rq3 <- draw_post(m_brms_error_hurdle_gamma,
"Error magnitude (pt)",
ERROR_COLOR,
PALETTE_ERROR,
adjust_in = 1,
breaks_x = seq(0, 10, by = 6))
g2_rq3 <- draw_post(m_brms_time_gamma,
"Response time (s)",
TIME_COLOR,
PALETTE_TIME,
adjust_in = 1,
breaks_x = seq(0, 18, by = 12))
g_combined_rq3 <- ggarrange(g1_rq3, g2_rq3, ncol = 2)
ggsave(g_combined_rq3, filename = "pic/RQ3-div.pdf", width = 6, height = 3.4, units = "in")
g_combined_rq3To investigate the effectiveness of each combination, we draw posterior samples from the model to get an estimate of all possible combinations. These posterior distributions show that how an average participant does with one combination, which corresponds to an experimental condition. Given our primary interest in error magnitude, we form a rank list of performance based on median of mean distribution.
global_order <- NULL
draw_all <- function( m,
title = "",
fill_color = "white" ,
palette_in,
adjust_in,
breaks_x = seq(0, 1, by = 0.5),
which = "error"
){
post_samples <- trials %>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
distinct() %>%
add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = FALSE) %>%
select(-c(".chain", ".iteration"))%>%
ungroup()%>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
mutate( # code it back so that we know which one is which...
Device = recode(Device, `0.5` = "HMD", `-0.5` = "Desktop"),
Motion = recode(Motion, `0.5` = "Dynamic", `-0.5` = "Static"),
Projection = recode(Projection, `0.5` = "Persp.", `-0.5` = "Orth."),
Shading = recode(Shading, `0` = "Simple", `0.5` = "A.O.", `-0.5` = "Flat"),
DataModel = recode(DataModel, `0.5` = "Image Data", `-0.5` = "Text Data"),
Dimensionality = recode(Dimensionality, `-0.5` = "2D Embedding", `0.5` = "3D Embedding"),
Condition = paste(Device, Motion, Projection, Shading, Dimensionality, DataModel, sep = seper))
b_CIs <- post_samples %>%
select(Condition, sample_value) %>%
group_by(Condition) %>%
median_qi() %>%
ungroup() %>%
mutate(lower95 = .lower, upper95 = .upper, median = sample_value,
s_lower95 = round(lower95, 2), s_upper95 = round(upper95, 2),
CIs = paste0(round(median, 2), ' [', s_lower95, ',', s_upper95, ']'))
if(which == "error"){
b_CIs <- b_CIs %>%
mutate(Condition = fct_reorder(Condition, .x = median, .desc = TRUE))
global_order <<- levels(b_CIs$Condition)
}else{
b_CIs <- b_CIs %>%
mutate(Condition = factor(Condition, levels = global_order))
}
median_max <- max(b_CIs$median)
median_min <- min(b_CIs$median)
g <- post_samples %>%
ggplot(aes(x = sample_value, y = Condition)) +
geom_vline(xintercept = median_max, linetype = "dashed", color = "black", size = LINE_SIZE) +
geom_vline(xintercept = median_min, linetype = "dashed", color = "black", size = LINE_SIZE) +
geom_richtext(data = b_CIs, mapping = aes(x = tail(breaks_x, n = 1) * 1.8,
y = Condition,
label = CIs),
label.color = NA,
color = TEXT_COLOR,
size = 3,
fill = NA,
hjust = 1,
label.padding = grid::unit(rep(0, 4), "pt")) +
geom_density_ridges(
fill = fill_color,
color = fill_color,
rel_min_height = 0.00005,
alpha = 0.5,
scale = .9,
size = 0.5
) +
scale_x_continuous(breaks = breaks_x, expand = expansion(mult = c(0.01, .8))) +
coord_cartesian(xlim = c(breaks_x[1], tail(breaks_x, n = 1))) +
ggtitle(title) + xlab("") + ylab("")
g
return(g)
}
g1_rq4 <- draw_all(m_brms_error_hurdle_gamma,
"Error magnitude (pt)",
ERROR_COLOR,
PALETTE_ERROR,
adjust_in = 1,
breaks_x = seq(0, 18, by = 6),
which = "error")
g2_rq4 <- draw_all(m_brms_time_gamma,
"Response time (s)",
TIME_COLOR,
PALETTE_TIME,
adjust_in = 1,
breaks_x = seq(0, 36, by = 12),
which = "time")
g_combined_rq4 <- ggarrange(g1_rq4, g2_rq4, ncol = 2)
ggsave(g_combined_rq4, filename = "pic/RQ4-all.pdf", width = 14, height = 15, units = "in")
ggsave(g1_rq4, filename = "pic/RQ4-error.pdf", width = 8, height = 12.5, units = "in")
ggsave(g2_rq4, filename = "pic/RQ4-time.pdf", width = 8, height = 12.5, units = "in")
g_combined_rq4This is the encoding to design the bubbles maps.
draw_tiles <- function( m,
title = "",
fill_color = "white" ,
stop_color = "white",
palette_in,
adjust_in,
breaks_x = seq(0, 1, by = 0.5),
which = "error"
){
post_samples <- trials %>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel) %>%
distinct() %>%
add_fitted_draws(m, value = "sample_value", scale = "response", re_formula = NA, allow_new_levels = FALSE) %>%
select(-c(".chain", ".iteration"))%>%
ungroup()%>%
select(Device, Motion, Projection, Shading, Dimensionality, DataModel, sample_value) %>%
mutate( Device = recode(Device, `0.5` = "2 HMD", `-0.5` = "1 Desktop"),
Motion = recode(Motion, `0.5` = "2 Dynamic", `-0.5` = "1 Static"),
Projection = recode(Projection, `-0.5` = "1 Orthographic", `0.5` = "2 Perspective"),
Shading = recode(Shading, `-0.5` = "1 Flat Shading", `0` = "2 Simple", `0.5` = "3 Ambient Occlusion"),
DataModel = recode(DataModel, `0.5` = "2 Image Data", `-0.5` = "1 Text Data"),
Dimensionality = recode(Dimensionality, `-0.5` = "2D Embedding", `0.5` = "3D Embedding"),
Condition = paste(Dimensionality, Device, Motion, DataModel, Projection, Shading, sep = seper))
b_CIs <- post_samples %>%
select(Condition, sample_value) %>%
group_by(Condition) %>%
median_qi() %>%
ungroup() %>%
mutate(lower95 = .lower, upper95 = .upper, median = sample_value,
s_lower95 = round(lower95, 3), s_upper95 = round(upper95, 3),
CIs = paste0(round(median, 3), ' [', s_lower95, ',', s_upper95, ']'),
index = row_number()
)
scale_factor <- 2.8
g <- b_CIs %>%
ggplot() +
theme(legend.position = "right", panel.grid.major = element_blank()) +
coord_fixed(ratio = 1) +
geom_tile(aes(x = 47, y = scale_factor * index), width = scale_factor, height = scale_factor, size = 0.1, fill = NA, color = "black") +
geom_circle(aes(color = median, x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(upper95)), fill = NA, size = 0.1) +
geom_circle(aes(fill = median, x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(median)), color = NA, size = 0.1) +
geom_circle(aes(x0 = 47, y0 = scale_factor * index, r = adjust_in * sqrt(lower95)), color = "white", fill = NA, size = 0.1) +
geom_text(aes(x = 20, y = scale_factor * index, label = Condition), size = 2) +
scale_x_continuous(limits = c(0, 50)) +
scale_fill_gradient(low = stop_color, high = fill_color, limits = c(breaks_x[1], tail(breaks_x, n = 1)), guide = "colourbar") +
scale_color_gradient(low = stop_color, high = fill_color, limits = c(breaks_x[1], tail(breaks_x, n = 1)), guide = "colourbar") +
scale_y_reverse() +
ggtitle(title) + xlab("") + ylab("")
g
return(g)
}
g1_rq4_tiles <- draw_tiles(m_brms_error_hurdle_gamma,
"Error magnitude (pt)",
ERROR_COLOR,
"white",
PALETTE_ERROR,
adjust_in = 0.35,
breaks_x = seq(3, 15, by = 4),
which = "error")
g2_rq4_tiles <- draw_tiles(m_brms_time_gamma,
"Response time (s)",
"#699923", #8bba45
"#f0fcd7",
PALETTE_TIME,
adjust_in = 0.24,
breaks_x = seq(0, 36, by = 12),
which = "time")
g_combined_rq4_tiles <- ggarrange(g1_rq4_tiles, g2_rq4_tiles, ncol = 2)
ggsave(g_combined_rq4_tiles, filename = "pic/RQ4-bubbles.pdf", width = 10.5, height = 19.5, units = "in")
g_combined_rq4_tilesdemographics <- metas %>% select(familiarML, familiarDR, familiarVIS, familiarVR)
demographics_stack <- stack(demographics)
demographics_stack$overall <- paste(demographics_stack$values, demographics_stack$ind, sep = " ")
demographics_table <- aggregate(data = demographics_stack, values ~ overall, length)
demographics_table$item <- sapply(demographics_table$overall, function(x){return(strsplit(x, split = " ")[[1]][2])})
demographics_table$option <- sapply(demographics_table$overall, function(x){return(strsplit(x, split = " ")[[1]][1])})
df <- data.frame(
item = c("normal", "glasses", "contacts", "others", "a little bit", "not at all"),
values = c(10, 5, 5, 2, 9, 23)
)
demo_color <- "#6959a6"
g_demographics <- ggplot() +
theme(legend.position = "top", panel.grid.major = element_blank()) +
#geom_tile(data = demographics_table, aes(x = option, y = item, fill = values)) +
geom_point(data = demographics_table, aes(color = values, x = as.numeric(option), y = item, size = sqrt(values)), stroke = 1) +
geom_point(data = df, aes(color = values, x = 8, y = item, size = sqrt(values)), stroke = 1) +
scale_fill_gradient(low = "#ffffff", high = demo_color, limits = c(0, 25), guide = "colourbar") +
scale_color_gradient(low = "#ffffff", high = demo_color, limits = c(0, 25), guide = "colourbar") +
scale_x_continuous(limit = c(0, 9), breaks = seq(1, 7, by = 1)) +
xlab("scale") + ylab("")
ggsave(g_demographics, filename = "pic/g_demographics.pdf", width = 3.5, height = 3.4, units = "in")
g_demographicsacc1 <- c(0.3978, 0.4054, 0.5998, 0.602, 0.605, 0.6052, 0.6528, 0.653, 0.6614, 0.6762, 0.6956, 0.7044, 0.713, 0.7162, 0.725, 0.73, 0.7326, 0.7362, 0.7398, 0.745, 0.7458, 0.7476, 0.751, 0.7538, 0.7668, 0.767, 0.7686, 0.7686, 0.769, 0.7744, 0.7748, 0.7754, 0.7758, 0.7778, 0.7788, 0.7794, 0.7806, 0.7828, 0.7834, 0.7848, 0.7858, 0.7882, 0.7888, 0.789, 0.7894, 0.7914, 0.7932, 0.7946, 0.7948, 0.7952, 0.7954, 0.7966, 0.7968, 0.7976, 0.799, 0.7992, 0.801, 0.8014, 0.8018, 0.8026, 0.8028, 0.8032, 0.804, 0.804, 0.8046, 0.805, 0.8056, 0.8058, 0.806, 0.8066, 0.8068, 0.8072, 0.8076, 0.808, 0.8084, 0.8098, 0.8108, 0.8108, 0.8108, 0.8118, 0.812, 0.8132, 0.8138, 0.8138, 0.8148, 0.816, 0.816, 0.8172, 0.8176, 0.8182, 0.8192, 0.8202, 0.8204, 0.8204, 0.8206, 0.8206, 0.8206, 0.8212, 0.8212, 0.822, 0.8224, 0.824, 0.8244, 0.8254, 0.8258, 0.8268, 0.8274, 0.8276, 0.8286, 0.8286, 0.829, 0.8302, 0.8314, 0.8314, 0.8326, 0.833, 0.8332, 0.8338, 0.834, 0.8352, 0.8358, 0.8364, 0.837, 0.8384, 0.8386, 0.8386, 0.8396, 0.8406, 0.8408, 0.8408, 0.8414, 0.842, 0.8428, 0.8436, 0.8442, 0.8444, 0.8444, 0.8446, 0.845, 0.8452, 0.847, 0.8472, 0.8472, 0.8474, 0.8474, 0.848, 0.85, 0.8504, 0.8506, 0.8514, 0.8516, 0.8518, 0.8534, 0.855, 0.8578, 0.858, 0.8582, 0.8584, 0.8584, 0.861, 0.8634, 0.8662, 0.8912, 0.8914, 0.8922, 0.8942, 0.8954, 0.896, 0.896, 0.896, 0.8962, 0.8964, 0.8966, 0.8968, 0.8972, 0.8974, 0.8978, 0.8978, 0.8978, 0.8982, 0.8982, 0.8984, 0.8984, 0.8986, 0.899, 0.8992, 0.8992, 0.8994, 0.8996, 0.8998, 0.8998, 0.8998, 0.8998, 0.8998, 0.9, 0.9, 0.9002, 0.9002, 0.9002, 0.9002, 0.9002, 0.9004, 0.9006, 0.9006, 0.9008, 0.901, 0.901, 0.9012, 0.9012, 0.9014, 0.9014, 0.9014, 0.9016, 0.9018, 0.9018, 0.9018, 0.902, 0.902, 0.902, 0.9022, 0.9024, 0.9026, 0.9026, 0.9026, 0.9026, 0.9026, 0.9028, 0.9028, 0.903, 0.903, 0.903, 0.903, 0.903, 0.9032, 0.9032, 0.9036, 0.9036, 0.9036, 0.9036, 0.9036, 0.9038, 0.9038, 0.9038, 0.9038, 0.9044, 0.9044, 0.9044, 0.9044, 0.9046, 0.9046, 0.905, 0.9052, 0.9052, 0.9052, 0.9054, 0.9054, 0.9056, 0.9056, 0.9056, 0.9058, 0.9058, 0.9058, 0.9058, 0.906, 0.906, 0.906, 0.906, 0.906, 0.9062, 0.9064, 0.9064, 0.9066, 0.9066, 0.9068, 0.9068, 0.9068, 0.9068, 0.907, 0.907, 0.9074, 0.9074, 0.9074, 0.908, 0.908, 0.9086, 0.9092, 0.9092, 0.9092, 0.9092, 0.9094, 0.9094, 0.9098, 0.9098, 0.91, 0.9106, 0.9114, 0.9114, 0.9126, 0.913, 0.9146)
acc2 <- c(0.192, 0.246, 0.394, 0.405, 0.427, 0.439, 0.476, 0.485, 0.488, 0.488, 0.49, 0.491, 0.491, 0.492, 0.492, 0.494, 0.496, 0.496, 0.497, 0.497, 0.497, 0.498, 0.498, 0.499, 0.499, 0.499, 0.5, 0.5, 0.5, 0.5, 0.5, 0.501, 0.504, 0.505, 0.505, 0.506, 0.506, 0.507, 0.507, 0.508, 0.508, 0.508, 0.508, 0.509, 0.509, 0.51, 0.51, 0.511, 0.511, 0.512, 0.512, 0.514, 0.515, 0.515, 0.515, 0.516, 0.522, 0.526, 0.53, 0.531, 0.531, 0.531, 0.534, 0.542, 0.548, 0.559, 0.564, 0.6, 0.625, 0.658, 0.681, 0.697, 0.722, 0.724, 0.726, 0.731, 0.735, 0.739, 0.739, 0.739, 0.739, 0.743, 0.747, 0.748, 0.753, 0.762, 0.764, 0.767, 0.776, 0.79, 0.79, 0.792, 0.795, 0.802, 0.803, 0.81, 0.813, 0.813, 0.817, 0.819, 0.82, 0.82, 0.823, 0.824, 0.824, 0.826, 0.828, 0.828, 0.834, 0.837, 0.838, 0.839, 0.841, 0.841, 0.842, 0.842, 0.843, 0.843, 0.844, 0.847, 0.848, 0.852, 0.853, 0.854, 0.854, 0.854, 0.856, 0.857, 0.857, 0.857, 0.859, 0.86, 0.862, 0.863, 0.863, 0.864, 0.865, 0.866, 0.871, 0.871, 0.874, 0.874, 0.875, 0.876, 0.879, 0.88, 0.881, 0.883, 0.883, 0.883, 0.884, 0.886, 0.886, 0.886, 0.887, 0.887, 0.888, 0.891, 0.891, 0.893, 0.893, 0.894, 0.896, 0.896, 0.896, 0.897, 0.898, 0.899, 0.901, 0.901, 0.901, 0.901, 0.902, 0.903, 0.904, 0.904, 0.904, 0.905, 0.905, 0.905, 0.907, 0.908, 0.908, 0.909, 0.91, 0.91, 0.91, 0.911, 0.911, 0.911, 0.911, 0.912, 0.912, 0.913, 0.913, 0.914, 0.917, 0.917, 0.92, 0.922)
acc3 <- c(0.838, 0.839, 0.8405, 0.8505, 0.8575, 0.8575, 0.865, 0.872, 0.872, 0.873, 0.8735, 0.8745, 0.882, 0.882, 0.8835, 0.885, 0.8865, 0.887, 0.887, 0.892, 0.892, 0.8935, 0.895, 0.8955, 0.896, 0.897, 0.8975, 0.898, 0.898, 0.8985, 0.8985, 0.899, 0.8995, 0.8995, 0.8995, 0.9, 0.9015, 0.902, 0.9025, 0.9035, 0.904, 0.9045, 0.9045, 0.905, 0.9055, 0.9055, 0.906, 0.906, 0.9075, 0.908, 0.9085, 0.9085, 0.9085, 0.9085, 0.9095, 0.9095, 0.9105, 0.9105, 0.9105, 0.911, 0.9115, 0.9115, 0.9115, 0.9115, 0.9125, 0.9125, 0.9125, 0.9125, 0.9125, 0.9125, 0.913, 0.9135, 0.9135, 0.9135, 0.914, 0.914, 0.914, 0.9145, 0.915, 0.915, 0.915, 0.9155, 0.9155, 0.9155, 0.9155, 0.916, 0.916, 0.916, 0.9165, 0.9165, 0.917, 0.917, 0.917, 0.917, 0.9175, 0.918, 0.918, 0.918, 0.918, 0.918, 0.918, 0.918, 0.9185, 0.9185, 0.9185, 0.9185, 0.9185, 0.919, 0.919, 0.919, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.9195, 0.92, 0.92, 0.92, 0.92, 0.9205, 0.9205, 0.9205, 0.9205, 0.921, 0.921, 0.9215, 0.9215, 0.9215, 0.9215, 0.9215, 0.922, 0.922, 0.922, 0.922, 0.922, 0.922, 0.9225, 0.9225, 0.9225, 0.9225, 0.923, 0.923, 0.923, 0.923, 0.9235, 0.9235, 0.9235, 0.9235, 0.9235, 0.9235, 0.924, 0.924, 0.924, 0.924, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.9245, 0.925, 0.925, 0.925, 0.925, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.9255, 0.926, 0.926, 0.926, 0.926, 0.926, 0.9265, 0.9265, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.927, 0.9275, 0.9275, 0.9275, 0.9275, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.928, 0.9285, 0.9285, 0.9285, 0.9285, 0.9285, 0.9285, 0.929, 0.929, 0.929, 0.929, 0.929, 0.929, 0.9295, 0.9295, 0.93, 0.93, 0.93, 0.93, 0.93, 0.9305, 0.9305, 0.9305, 0.931, 0.9315, 0.9315, 0.9315, 0.9315, 0.932, 0.9325, 0.9325, 0.933, 0.9335, 0.9335, 0.9335, 0.934, 0.9345, 0.9345, 0.936, 0.936, 0.937, 0.937, 0.94, 0.9405)
bars_color <- "#6959a6"
g_bar1 <- ggplot() +
theme(panel.grid.major.x = element_blank()) +
geom_linerange(aes(ymax = 1, ymin = 0, x = acc1), size = 0.07, color = bars_color) +
scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0)) +
ggtitle("CIFAR-10")
g_bar2 <- ggplot() +
theme(panel.grid.major.x = element_blank()) +
geom_linerange(aes(ymax = 1, ymin = 0, x = acc2), size = 0.07, color = bars_color) +
scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0))+
ggtitle("bABI")
g_bar3 <- ggplot() +
theme(panel.grid.major.x = element_blank()) +
geom_linerange(aes(ymax = 1, ymin = 0, x = acc3), size = 0.07, color = bars_color) +
scale_x_continuous(limits=c(0,1), breaks = seq(0, 1, by = 0.1), expand=c(0,0)) +
scale_y_continuous(limits=c(0,1), breaks = seq(0, 1, by = 1), expand=c(0,0)) +
ggtitle("F-MNIST")
g_bars <- ggarrange(g_bar1, g_bar2, g_bar3, ncol = 1, align="v")
ggsave(file = "pic/acc_bar.pdf", g_bars , width = 4, height = 2, units = "in")
g_bars